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ABSTRACT 

A hybrid  turbulence  modeling  approach 
combining  unsteady  Reynolds-averaged  Navier- 
Stokes  (URANS)  equation-based  modeling  and  large 
eddy  simulation  (LES)  is  applied  to  turbulent  free- 
surface  flows  around  surface-piercing  bodies 
involving  wave-breaking  and  bubble  formation. 
Based  on  a two-equation  k-s  turbulence  model,  the 
hybrid  model  reduces  to  a RANS  turbulence  model 
or  a subgrid-scale  turbulence  model  depending  on  the 
local  mesh  resolution  (filter  size)  and  the  local 
integral  turbulence  length-scale.  The  computations 
were  carried  out  using  a second-order-accurate  finite- 
volume  Navier-Stokes  solver  that  permits  use  of 
arbitrary  unstructured  meshes  and  local  mesh 
refinement.  Volume-of-fluid  (VOF)  method  was 
employed  to  capture  the  free-surface.  The  URANS- 
LES  hybrid  approach  is  shown  to  reproduce  the 
salient  features  of  the  turbulent  free-surface  flows 
around  the  surface-piercing  bodies  with  reasonable 
accuracy  on  relatively  coarse  meshes  on  which  the 
fidelity  of  LES  degrades. 

INTRODUCTION 

One  of  the  difficulties  encoimtered  in 
numerically  predicting  turbulent  free-surface  flows 
comes  from  turbulence  modeling.  Turbulent  structure 
near  free-surface  is  markedly  different  from  those  in 
regions  away  from  free-surface.  Anisotropy  of 
turbulence  and  energy  backscatter  are  among  much 
discussed  features  of  free-surface  turbulence  (e.g., 
Shen  and  Yue,  2001).  In  this  respect,  free-surface 
turbulence  shares  its  characteristics  with  near-wall 
turbulence.  And  yet,  the  physics  of  free-surface 
turbulence  widely  varies  depending  on  Froude  number 
and  Weber  number  (see  Brocchini  and  Peregrine, 
2001).  For  violent  free-surface  flows  around  surface- 
piercing bodies  involving  wave-induced  flow 
separation,  wave-breaking,  and  bubble  formation,  it  is 


a formidable  task  to  gain  even  a qualitative 
understanding  of  interactions  among  surface  waves, 
boundary  layer,  wake,  and  their  effects  on  free-surface 
turbulence,  let  alone  to  find  any  models  reflecting  those 
effects.  Adds  to  all  this  the  difficulty  of  modeling  the 
physics  of  turbulence  associated  with  interfacial 
transfer  of  mass,  momentum,  and  kinetic  energy 
occurring  in  two-phase  flows.  Given  the  challenges,  it 
is  not  difficult  to  see  why  Reynolds-averaged  Navier- 
Stokes  (RANS)  equations-based  turbulence  modeling 
is  considered  less  than  adequate  for  this  class  of  free- 
surface  flows. 

Large  eddy  simulation  (LES)  has  been 
considered  to  be  fundamentally  more  suitable  for  the 
subject  flow,  in  view  of  the  dominant  role  played  by 
large-scale  turbulent  structures.  Dommermuth  and 
Novikov  (1993),  Shen  and  Yue  (2001),  Kawamura  el 
al.  (2002),  Yue  el  al.  (2005),  and  most  recently  Kim 
and  Rhee  (2006)  are  among  those  who  attempted  LES 
for  simulation  of  turbulent  free-surface  flows.  For  all 
the  merits  of  LES  for  flows,  LES  is  computationally 
intensive,  often  prohibitively.  Despite  the  remarkable 
increase  in  computational  power  in  terms  of  hardware 
and  software.  Alleviating  the  cost  of  LES  is  an  area  of 
active  research  these  days.  The  first  fruit  of  such 
efforts  came  in  the  so-called  Detached  Eddy 
Simulation  (DES)  approach  proposed  by  Spalart  el  al. 
(1997).  hi  broad  terms,  DES  sets  out  to  selectively  use 
either  RANS  or  subgrid-scale  (SGS)  turbulence 
modeling  in  the  flow  domain  in  question  based  on 
certain  criterion.  The  turbulence  model  adopted  in 
Spalart’s  DES  approach  reduces  to  the  Spalart  and 
Allmaras’  eddy-viscosity  transport  (RANS)  model  in 
near-wall  region,  whereas  it  reduces  to  a subgrid-scale 
(SGS)  turbulence  model  away  from  wall.  This  earliest 
version  of  DES  approaches  allows  one  to  use  highly 
stretched  meshes  in  attached  boundary  for  use  with  the 
RANS  modeling,  which  significantly  saves  cell  counts 
and,  as  a consequence,  reduce  the  computational  cost 

The  present  study  is  a sequel  to  the  earlier 
LES  study  by  the  authors  (Kim  and  Rhee,  2006). 
What  we  tried  to  do  here  in  this  study  is  to  critically 


evaluate  a URANS/LES  hybrid  turbulence  modeling 
approach  for  turbulent  free-surface  flows  around 
surface-piercing  bodies.  The  hybrid  model  is  based  on 
two-equation  k-e  turbulence  model.  The  same 
numerical  methods  and  algorithms  were  used  as  in  the 
previous  study.  Free-surface  is  resolved  with  volume- 
of-fluid  (VOF)  method.  The  same  surface-piercing 
hydrofoil  is  considered  with  a rectangular  plan-form 
and  NACA-0024  cross-section  experimentally  studied 
in  the  towing  tank  at  the  University  of  Iowa  (Zhang 
and  Stem,  1996;  Metcalf  et  al.,  2006).  We  considered 
the  case  of  Fr  = 0.37  in  this  study.  The  corresponding 
Reynolds  numbers  (Re  = U0c/v)  are  1.52  x 106.  In 
addition,  the  flow  around  a surface-piercing  circular 
cylinder  was  newly  considered  in  this  study,  which 
was  experimentally  studied  by  Inoue  et  al.  (1993). 
LES  as  well  as  DES  computations  have  been  carried 
out  for  the  circular  cylinder  case  at  the  Froude  number 
(/•>)  of  0.8,  with  the  corresponding  Reynolds  number 
ofTte  = 2.7  x 104. 

The  computations  have  been  carried  out  using 
a parallelized,  finite-volume  Navier-Stokes  solver  that 
is  second-order  accurate  in  time  and  space. 


GOVERNING  EQUATIONS) 


The  governing  equations  solved  in  DES  are 
a single  set  of  the  filtered  (averaged)  Navier-Stokes 
equations  that  are  shared  by  gas  (air)  and  liquid 
(water)  phases  or  their  mixture. 

The  filtered  Navier-Stokes  equations  can  be 
written  as 


and 
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where  the  overbar  denotes  the  filtering  operation,  p 
the  piezometric  pressure  ( p = ps+  pgZ ). 

The  presence  of  two  phases  (air  and  water) 
in  the  computational  domain  is  represented  by 
volume  fraction,  ;<x,t)  (Harlow  and  Welch,  1965). 
Volume-fraction  in  intermittent  two-phase  flows  can 
be  interpreted  as  a volume-averaged  intermittency 
function,  I(x,t)  whose  value  takes  either  0 or  1 
depending  on  whether  the  point  (x)  in  question  is  in 
water  or  air.  Thus, 


y(x, /)  = /(x,  t)  = y jj£  /(x,  /)  dV  (3) 

V c c 

Note  that  the  volume-fraction  changes  in  both  space 
and  time,  whose  value  ranges  between  0 and  1. 
Volume-fraction  is  obtained  from  its  transport 
equation 


dy 

~8t 


(4) 


The  physical  properties  such  as  density  and 
dynamic  viscosity  are  computed  as  functions  of 
/from 


p = (l-y)pw+ypa  (5) 

p = {l-y)pw  +ypa 

where  the  subscripts  “a”  and  “w”  denote  air  and 
water,  respectively. 

Insofar  as  the  fluid  density  fluctuates  in  the 
vicinity  of  free-surface  with  the  volume-fraction,  the 
filtered  variables  in  Equation  (1)  and  Equation  (2), 
except  the  density,  can  be  interpreted  as  density- 
weighted  (Favre-averaged)  variables  as  employed  for 
compressible  flows. 

NUMERICAL  METHOD 

The  computations  were  carried  out  using  a 
development  version  of  FLUENT  V6.3.  FLUENT 
employs  a cell-centered  finite-volume  method  based 
on  a multi-dimensional  linear  reconstruction  scheme 
that  permits  use  of  computational  elements  (cells)  with 
arbitrary  polyhedral  cell  topology  including 
quadrilateral,  hexahedral,  triangular,  tetrahedral, 
pyramidal,  prismatic,  and  hybrid  meshes.  The 
solution  gradients  at  cell  centers  are  obtained  by 
applying  Green-Gauss  theorem  based  on  a node-based 
quadrature  (Kim  et  al.,  2003).  In  the  present  study, 
the  discretization  schemes  and  the  solution  algorithms 
were  judiciously  chosen  to  ensure  accuracy,  efficiency 
and  stability  of  the  numerical  solutions,  all  of  which 
are  essential  for  successful  LES.  For  discretization  of 
convection  terms  in  the  filtered  momentum  equations, 
we  adopted  a bounded  central-differencing  (BCD) 
scheme  that  is  essentially  non-dissipative,  second- 
order-accurate  central-differencing  scheme  combined 
with  a nonlinear  flux-limiter  activated  only  when  and 
where  needed  to  suppress  numerical  oscillations.  For 
discretization  of  advection  term  in  the  VOF  equation, 
we  employed  a modified  high  resolution  interface- 
capturing (HRIC)  scheme  that  is  compressive  and 
therefore  ideally  suited  for  resolving  sharp  interface. 
Surface  tension  effect  is  included  using  the  continuum 


surface  force  (CSF)  model  that  effectively  represents 
surface  tension  as  a volumetric  source  (forcing)  term 
in  the  momentum  equations. 

The  system  of  discretized  governing 
equations  is  solved  using  point-implicit  Gauss-Seidel 
relaxation  along  with  algebraic  multi-grid  (AMG) 
method  to  accelerate  solution  convergence.  The  N-S 
solver  and  the  SGS  turbulence  model  arc  hilly 
parallelized  using  MPI. 

An  implicit  fractional-step  method  (Kim  and 
Makarov,  2005)  in  combination  with  a second-order 
accurate,  three-level  backward-differencing  scheme 
for  time-discretization  was  employed  to  advance  the 
solution  in  time.  Unlike  explicit  schemes,  implicit 
time-advancement  scheme  allows  one  to  take  as  large 
a time-step  size  as  needed  to  guarantee  a desired 
temporal  resolution,  bringing  a considerable  speed-up 
in  simulations  of  transient  flows.  In  this  algorithm,  the 
momentum  equations  are  decoupled  from  the 
continuity  equation  using  an  approximate  factorization 
of  the  coupled  Navier-Stokes  equations. 

In  the  present  computations,  the  volume- 
fraction  equation  was  solved  first  in  the  beginning  of 
each  time-step,  which  is  followed  by  the  remaining 
steps  with  the  newly  updated  physical  properties  of 
the  mixtiue. 

TURBULENCE  MODELING 

The  DES  method  originally  proposed  by 
Spalart  et  al.  (1997)  has  been  widely  used  for 
aerodynamics  applications  (Nikitin  et  al.,  2000; 
Travin  et  al.,  2000;  Cokljat  and  Liu,  2002;  Spalart  et 
al.,  2003).  An  interesting  study  had  been  conducted 
for  an  A-airfoil  for  the  LESFOIL  European  project. 
The  DES  study  of  Cokljat  and  Liu  (2002)  showed  an 
excellent  agreement  despite  the  relatively  coarse 
mesh  used  for  the  DES  computation.  To  obtain  a 
similar  acciuacy  on  the  same  flow  using  LES,  a 
substantially  larger  mesh  had  to  be  used  (Mary  and 
Sagaut,  2002). 

In  spite  of  its  encouraging  performance  for 
massively  separated  flows,  DES  for  attached  flows, 
e.g.  channel  flows  (e.g.  Nikitin  et  al.,  2000),  showed 
a fairly  strong  dependency  on  grid  resolution.  In  that 
work,  the  skin-friction  coefficient  was  underpredicted 
by  about  15  %.  Furthermore,  DES  simulations  on 
separated  flows  as  studied  by  Travin  et  al.  (2000) 
also  revealed  various  short  comings.  One  of  the  main 
drawbacks  of  DES  approaches  in  general  lies  with 
the  interface  region  between  RANS  zone  and  LES 
zone  where  two  completely  different  turbulence 
models  are  to  be  matched. 

Strelets  (2000)  performed  DES 
computations  of  several  massively  separated  flows 
using  a new  hybrid  model  based  on  the  k-oi  SST 


model.  This  two-equation-based  hybrid  model 
offers  two  potential  advantages  compared  to  the 
original  DES  approach  of  Spalart  et  al.  (1997). 
Firstly,  modem  two-equation-based  RANS  models 
provide  better  accuracy  than  the  one-equation  model 
- or  lower-order  turbulence  models  - for  a wide  range 
of  flows  including  boundary  layers,  wakes,  and  other 
free  shear  flows.  Thus,  the  DES  approach  based  on 
advanced  two-equation  turbulence  models  would 
benefit  from  the  superior  performance  of  the 
advanced  RANS  models.  Secondly,  unlike  the  one- 
equation  based  DES  model,  the  demarcation  of 
RANS  and  LES  regions  in  the  two-equation-based 
DES  approach  is  not  any  more  based  on  the 
proximity  to  wall  (wall  distance)  only.  Two-equation 
RANS  turbulence  models  such  as  k-e  or  k-co  models 
give  an  estimate  of  integral  length-scale  of  turbulence 
that  can  be  used  to  determine  whether  to  invoke 
subgrid-scale  turbulence  model  or  RANS  turbulence 
model.  The  decision  is  based  on  upon  the  integral 
length-scale  and  the  grid-filter  size. 

Among  many  hybrid  RANS/LES  models 
available,  we  picked  the  DES  approach  based  on 
realizable  k-e  turbulence  model  Shih  et  al.  (1995). 
The  modeled  transport  equation  for  k in  this  model  is 
given  by: 


where  Pk  is  the  production  of  turbulent  kinetic 
energy,  and  the  last  term  on  r.h.s.  of  the  above 
equation  is  modified  dissipation  rate  formulation 
used  in  conjunction  with  the  present  DES  approach. 
The  turbulent  length  scale  is  defined  as: 

tDES  =MIN(k3l2/S,CDES  A)  (7) 

where  A is  the  grid  length-scale.  CDES  is  calibration 

constant  and  the  value  of  0.61  is  adopted  in  the 
present  study  and  A is  the  grid  length  scale 
computed  from  a = k[V  where  V is  the  volume  of  the 
cell.  The  transport  equation  for  e was  taken  from  the 
realizable  k-e  equation  of  Shih  et  al.  (1995). 

For  the  LES  computations  for  the  circular 
cylinder  case,  the  same  dynamic  Smagorinsky  model 
was  used  as  the  one  employed  in  the  earlier  study  of 
the  authors  (Kim  and  Rhee,  2006). 

COMPUTATIONAL  DOMAIN  AND  MESH 


Circular  cylinder 

The  upstream  inlet  boundary  and  the 
downstream  exit  boundary  were  placed  at  8.5/J  and 
20.5 D forward  and  aft  from  the  cylinder  center, 
respectively.  The  lateral  boundary  is  5.0/J  away 
from  the  cylinder  axis.  The  top  and  the  bottom 
boundaries  of  the  domain  were  placed  at  2.0/7  and 
5.0D  above  and  below  the  undisturbed  free-surface, 
respectively.  The  blockage  ratio  of  the  cylinder 
(H/D,  where  H is  the  height  of  the  domain)  is 
approximately  4.8  %. 

In  order  to  quantify  the  mesh-dependency,  the 
computations  were  made  with  two  meshes  with  1.5 
million  and  2.7  million  elements. 

Hydrofoil 

The  computational  domain  is  bounded  by  an 
upstream  inlet,  a downstream  exit,  top  and  bottom 
planes,  and  a pair  of  far-field  boundaries  in  the  lateral 
direction.  The  upstream  inlet  boundary  and  the 
downstream  exit  boundary  were  placed  at  3.0  chord- 
lengths  (c)  and  8.54  chord-lengths  forward  and  aft 
from  the  leading-edge  of  the  hydrofoil,  respectively. 
The  top  and  the  bottom  boimdaries  of  the  domain 
were  placed  at  0.21  c and  2.5  c above  and  below  the 
undisturbed  free-surface,  respectively.  The  lateral 
boundaries  are  at  5.0  c from  the  foil  center-plane. 

In  view  of  the  relatively  simple  geometry,  block- 
structured  hexahedral  meshes  with  a C-H  topology 
were  used  in  this  study.  Two  meshes  were  used  in 
this  study.  The  coarse  mesh  has  1.8  million 
elements,  and  the  fine  mesh  has  4.7  million  cells.  The 
fine  mesh  was  generated  by  locally  refining  all  the 
cells  in  the  coarse  mesh  in  the  strip  between  z = - 
0.15c  and  z = 0.1c.  The  near-wall  mesh  resolution 
for  the  4.7  million  cell  mesh  is  such  that  the  distance 
from  the  cylinder  surface  at  the  wall-adjacent  cells  is 
in  the  neighborhood  of  3 x 10"4  chord-length,  which 
translates  to  y+  in  the  range  of  3 ~ 35,  yet  staying 
close  to  y+  = 10  on  most  of  the  foil  surface.  The  fine 
mesh  has  twice  finer  than  the  coarser  mesh  in  the 
near-wall  region. 

BOUNDARY  CONDITIONS  AND  OTHER 
COMPUTATIONAL  DETAILS 

For  both  the  circular  cylinder  and  the 
hydrofoil  case,  on  the  upstream  and  lateral  inlet 
boundaries,  total  pressure  is  specified.  On  these 
“pressure-inlet”  boundaries  supported  in  FLUENT, 
velocity  is  determined  from  the  specified  total 
pressure  using  the  dynamic  pressure  and  the 
hydrostatic  pressure  at  given  vertical  location,  and 
volume  fraction  values  are  simply  extrapolated.  At 
the  downstream  exit  boundary,  the  hydrostatic 


pressure  is  determined  from  the  depth  below  the  free- 
surface,  while  other  variables  including  volume- 
fraction  are  extrapolated.  On  the  top  and  bottom 
boundaries,  symmetry  condition  is  imposed.  On  the 
foil  surface,  we  employed  a generalized  wall-fimction 
approach  that  invokes  proper  wall-laws  depending  on 
the  y+  value  at  wall-adjacent  cells.  Thus,  the  wall 
model  is  applicable  to  the  entire  inner  layer  including 
the  viscous  sub-layer,  buffer  zone,  and  logarithmic 
layer. 

The  non-dimensional  time-step  size  (At) 
used  for  the  computations  is  0.0002  sec.  The  time- 
step  size  was  determined  based  on  the  estimates  of 
the  characteristic  length  and  time  scales  of  the  eddy 
size  to  be  resolved.  The  smallest  time-scale  to  be 
resolved  was  computed  from r = £/u' , where  l was 
taken  as  0.02  c,  and  u'  as  0.5 U0.  From  these  rough 
estimates,  one  turnover  of  the  smallest  resolved 
eddies  will  be  resolved  with  roughly  200  time  steps. 

To  shorten  the  initial  transient  period  of  the 
solution  and  to  quickly  attain  statistically  stationary 
state,  we  took  steady  RANS  solution  as  the  initial 
condition. 

RESULTS 

Circular  cylinder 

The  results  for  the  circular  cylinder  case  obtained 
with  the  fine  mesh  (2.7M  cell)  mesh  for  Re  = 2.8  x 
104,  and  Fr  = 0.8  are  discussed  here.  An  overall 
impression  of  the  free-surface  around  the  circular 
cylinder  predicted  by  the  DES  is  shown  in  Figure  1, 
along  with  that  obtained  using  LES  reported  earlier 
(Kim  and  Rhee,  2006).  The  bow  wave  and  the 
downstream  wave  pattern  arising  from  the  interaction 
between  the  turbulent  wake  and  the  free-surface  are 
well  captured  by  the  DES.  As  evident  in  the  figures, 
a vigorous  aeration  (ventilation)  takes  place  in  the 
near- wake.  As  observed  in  the  LES,  large  air  bubbles 
were  seen  to  be  entrained  down  to  nearly  0.1D  below 
the  undisturbed  free-surface  in  the  DES  results.  Some 
of  the  vertically  oriented  vortices  leave  whirls  on  the 
free-surface  with  conspicuous  depression  around 
their  cores.  Overall,  the  DES  results  closely 
resemble  the  LES  results.  However,  the  flow 
structures  predicted  by  DES  appear  to  be  slightly 
more  coarse-grained  than  what  the  LES  yields.  In 
addition,  the  depression  of  the  free-surface  predicted 
by  the  DES  seems  smaller  than  the  LES  prediction. 

Figure  2 shows  the  contours  of  time-averaged 
elevation  of  the  free-surface  (y  = 0.5)  computed  from 
the  experiment  (Inoue  el  al.,  1993).  In  Figure  3 are 
depicted  the  time-averaged  free-surface  elevations 
predicted  by  the  DES  and  the  LES.  It  is  seen  that  the 


main  features  of  the  free-surface  elevation  are  well 
captured  by  the  DES  computation  such  as  the  wave 
pattern  and  height  in  front  of  the  cylinder  and  in  the 
wake.  However,  the  contours  clearly  show  that  the 
depression  of  the  free-surface  predicted  by  DES  is 
appreciably  smaller  than  the  LES  prediction, 
especially  in  the  near-wake  region  and  in  the 
shoulder  wave  region. 


Figure  4 are  shows  the  time-averaged  wave 
profiles  in  the  spanwise  (y)  direction  at  two  axial 
locations  (x/D  = 0.9,  2.0),  predicted  by  the  DES  and 
the  LES.  Considering  that  the  free-surface  violently 
fluctuate,  and  as  a result,  the  width  of  the  surface- 
layer  becomes  large,  the  predicted  profiles  are  plotted 
for  the  volume-fraction  values  of  f = 0.3,  0.4,  0.5. 
The  large  band-width  of  the  time-averaged  surface 
elevation  for  the  range  of  volume-fraction,  especially 
in  the  near-wake  of  the  cylinder  ( x/D  < 1.0),  is  an 
indication  of  large  fluctuation  of  the  free-surface  and 
vigorous  ventilation  occurring  near  the  free-surface. 
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Figure  3.  Contours  of  the  time-averaged  free-surface 
elevation  (Re  = 2.8  x 104,  Fr  = 0.8):  (a)  DES 
prediction;  (b)  LES  prediction 


Figure  1.  Instantaneous  free-surface  generated  by 
the  flow  around  a surface-piercing  circular  cylinder, 
colored  by  the  surface  elevation  (Re  = 2.8  x 104,  Fr  = 
0.8) 
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Figure  4.  Time-averaged  free-surface  elevation  (Re 
= 2.8  x 104,  Fr  = 0.8)  at  two  axial  locations  (x/D  = 
0.9,  2.0).  Note  that  the  predictions  for  three  different 
values  of  volume-fraction  are  shown. 

Figure  4 show  that  the  DES  yields  substantially 
smaller  depression  of  the  free-surface  in  the  wake 
region  than  the  LES  prediction.  One  observation 
worthy  of  mentioning  here  with  regard  to  the  smaller 
depression  of  the  free-surface  in  the  near-wake 
predicted  by  the  DES  is  that  the  DES  prediction  also 
gave  a depth-averaged  drag  coefficient  of  around  0.5, 
a much  smaller  value  than  the  LES  prediction  (CD  = 
0.95)  and  the  experimental  correlation  ( CD  = 1.0). 
Although  not  shown  here,  the  DES  prediction  yielded 
a substantially  narrower  wake  than  the  LES 
prediction.  In  view  of  the  fact  that  DES  model 
reduces  to  RANS  turbulence  model  in  the  near-wall 
region,  the  RANS  turbulence  model  (realizable  k-s) 
used  in  this  particular  DES  approach  is  most  likely  to 
blame  for  the  discrepancy  observed.  Indeed, 
efficacy  of  DES  - URAS/LES  hybrid  approaches  - 
would  surely  hinge  upon  the  performance  of  RANS 
turbulence  models  they  adopt. 


(b)  LES 


Figure  5.  Contours  of  the  r.m.s.  vertical  velocity 
fluctuation  on  the  free  surface  (Re  = 2.8  x 104,  Fr  = 
0.8):  top  - predicted  by  the  present  DES  and  LES 


Figure  5 depicts  the  contours  of  the  r.m.s. 
vertical  (z-)  velocity  fluctuation  plotted  on  the  mean 
free-surface  predicted  by  the  DES  and  the  LES 
computations,  which  give  an  useful  information  on 
the  rums,  free-surface  elevation  fluctuation.  What 
these  contour  plots  indicate  is  that  the  DES 
considerably  dampens  the  normal  velocity  fluctuation 
in  the  wake  in  comparison  to  the  LES  prediction.  As 
we  will  see  later,  this  is  also  the  case  with  the 
hydrofoil  case. 

In  summary,  the  DES  appears  to  predict  the 
gross  features  of  the  turbulent  free-surface  flow 
around  a surface-piercing  circular  cylinder  with  a 
reasonable  accuracy.  However,  it  tends  to  under- 
predict the  depression  of  the  free-surface  in  the  wake 
and  the  level  of  the  fluctuation  of  the  free-surface. 


(b)  LES 

Figure  6.  Instantaneous  surface  elevation  around  the 
surface-piercing  hydrofoil  predicted  by  the 
RANS/LES  hybrid  model  and  LES  at  Fr  = 0.37 


An  overall  impression  of  the  free-surface 
around  the  surface-piercing  hydrofoil  predicted  using 
the  present  DES  approach  is  shown  in  Figure  6 along 
with  the  LES  prediction  for  comparison.  The 
instantaneous  surface  elevations  were  visualized 
using  iso-surface  of  volume-fraction  for  y = 0.5.  The 
main  characteristics  of  the  free-surface  portrayed  in 
these  figures  closely  resemble  what  is  shown  in  the 
photographs  taken  in  the  experiment  (Figure  6 in 
Metcalf  el  al.,  2006).  They  are  also  qualitatively 
similar  to  those  observed  by  Pogozelski  and  Katz 
(1997)  with  the  free-surface  around  a surface- 
piercing strut.  The  main  difference  between  the  DES 
and  the  LES  predictions  is  the  degree  of  details  of  the 
flow  resolved  by  the  predictions.  The  LES  prediction 
appears  to  capture  finer  scales  of  motion  compared  to 
the  DES  result,  most  likely  because  the  DES  model 


(b)  LES 

Figure  7.  Ratio  of  turbulent  iscosity  to  molecular 
viscosity  (vt/v)  - plotted  in  the  range  of  0 < 
vt/v  <1000 


In  Figure  7 is  shown  the  comparison  of  the  ratio  of 
turbulent  viscosity  to  molecular  viscosity  at  an 
instance  of  time.  Not  surprisingly,  the  DES  approach 
yields  much  higher  - by  up  to  roughly  two  orders  of 
magnitude  - turbulent  viscosity  in  regions  away  from 
the  hydrofoil.  This  is  an  indication  that  the  DES 
effectively  reverts  back  to  RANS  turbulence  model 
(realizable  k-e  model  in  this  study)  in  some  parts  of 
the  flow  domain. 

How  the  solution  domain,  at  an  instance  of 
time,  is  demarcated  into  RANS  and  LES  regions  is 
visualized  in  Figure  8 using  a color  map  constructed 
using  the  size  of  the  integral  length  scale  of 
turbulence  relative  to  the  length-scale  of  the  grid.  In 
the  figure,  the  RANS  region  is  shown  in  blue. 


whereas  the  LES  region  is  depicted  in  red.  With  a 
caveat  in  mind  that  the  demarcation  is  shown  near  the 
free-surface,  the  color  map  indicates  that  the  RANS 
region  resides  in  the  thin  boundary  layer  region  near 
the  hydrofoil  and  in  the  region  far  way  from  the 
hydrofoil.  The  LES  region  is  between  the  two 
regions.  One  important  point  to  stress  here  is  that  the 
RANS  region  in  the  far  field  could  have  not  been 
recovered  by  the  original  DES  of  Spalart  el  al. 
(1997)  which  is  based  on  the  relative  size  of  wall 
distance  to  grid  spacing  as  the  only  demarcation 
criterion. 


Figure  8.  Color  map  showing  demarcation  of  RANS 
and  LES  region  - The  RANS  region  is  shown  in  blue, 
and  the  LES  region  in  red. 


Figure  9 shows  the  contours  of 
instantaneous  longitudinal  vorticity  component  at  a 
cross-flow  plane  (x/c  = 0.7)  predicted  by  the  DES  and 
the  LES.  The  free-surface  is  also  shown  in  these 
figures  in  solid  lines  near  the  top  of  each  figiure. 
These  contour  plots  shed  light  on  the  strength  and  the 
length  scale  of  the  subsurface  vortical  structure.  The 
LES  (Kim  and  Rhee,  2006),  shown  in  Figure  9(b) 
reproduce  the  characteristic  features  of  the  subsurface 
vortical  structures  such  as  the  extent  and  the  length- 
scales  of  the  vortices  and  the  occurrence  of  pairs  of 
counter-rotating  vortices  observed  by  Pogozelski  and 
Katz  (1997).  The  DES  prediction  appears  to  be 
qualitatively  similar  to  the  vortical  structures  pictured 
by  the  LES.  However,  the  LES  prediction  exhibits  a 
far  richer  feature  than  the  DES  result.  In  addition,  the 
vorticity  predicted  by  the  DES  is  much  weaker  than 
the  LES  prediction. 
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(b)  LES  (Kim  and  Rhee,  2006) 


Figure  9.  Contour  of  instantaneous  axial  vorticity 
component  ( cox)  between  -100  <a>x  < 100  viewed 
from  upstream 


Figures  10  show  the  contours  of  the  x- 
component  of  the  time-averaged  wall-shear  stress  on 
the  hydrofoil  surface  predicted  by  the  DES,  which 
clearly  indicates  that  there  is  a pocket  of  flow 
reversal  (liilqlcr  the  mean  free-surface  - a wedge- 
shaped  region  below  the  mean  free-surface,  colored 
in  blue  in  the  figure.  The  LES  reported  earlier  also 
captured  the  wedge-shaped  region  of  flow  reversal. 
The  contour  of  the  wall-shear  indicates  that  the  flow 
underneath  the  pocket  of  separation,  being 
accelerated  past  the  bow  and  the  shoulder  driven  by 
strong  favorable  pressure  gradient,  penetrates  deeply 
under  the  free-surface  and  remains  attached  nearly  to 
the  trailing-edge.  The  extent  and  the  topological 
structure  of  the  flow  separation  are  largely  consistent 
with  those  described  by  Zhang  and  Stem  (1996)  and 
Pogozelski  and  Katz  (1997). 


Figure  10.  Contour  of  the  time-averaged  x- 
component  of  the  wall-shear  stress  on  the  hydrofoil 
surface  predicted  by  the  present  DES 


Figure  1 1 shows  the  contour  of  the  time- 
averaged  free-surface  elevation  predicted  by  the  DES 
on  the  fine  mesh  (4.7  million  cell),  along  with  the 
measurement  (Metcalf  et  al.  2006).  The  DES 
prediction  closely  matches  the  measured  elevation. 
The  Kelvin  waves  away  from  the  hydrofoil  are 
closely  reproduced  by  the  DES  prediction,  in  terms  of 
the  height  and  the  position  of  the  wave  crests  and 
troughs  and  the  divergence  angle  of  the  waves.  The 
DES  apparently  overpredicts  the  surface  elevation  in 
the  central  portion  of  the  near-wake,  as  compared  to 
the  measurement.  Nonetheless,  the  present  DES 
result  is  quite  promising. 
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Figure  11.  Time-averaged  free-surface  wave 
elevation  predicted  by  the  present  DES  at  Fr  = 0.37 
Top  - prediction  (fine  mesh);  Bottom  - measured. 
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Figure  12.  Time-averaged  wave  profiles  along  the 
hydrofoil  foil  surface  predicted  by  (a)  DES;  (b)  LES 


Figures  12(a)  and  12(b)  depict  the  wave 
profiles  along  the  hydrofoil  surface  predicted  by  the 
DES  and  the  LES,  respectively,  along  with  the  ones 
determined  experimentally  (Zhang  and  Stem,  1996; 
Metcalf  et  al.,  2006).  The  time-averaged  wave 
profiles  are  shown  in  Figures  for  a few  different 
values  of  volume-fraction;  y = 0.2,  0.3,  0.4,  and  0.5. 
It  should  be  stressed  that  the  large  band-width  of  the 
wave  elevations  for  the  mean  volume-fraction 
ranging  from  y = 0.2  to  y = 0.5  is  an  evidence  of  the 
large  fluctuation  of  wave  elevation  and  the 
occurrence  of  intermittent  two-phase  flow  in  the 
separated  region  and  in  the  near-wake.  It  is  seen  that, 
in  terms  of  the  time-averaged  wave  profiles  along  the 
body  surface  (near-field),  the  DES  results  closely 
match  the  LES  predictions.  The  mean  wave  profile 
measured  by  the  servo-driven  wave  gage  differs 
significantly  from  those  determined  by  the  visual 
observations.  The  DES  and  the  LES  predictions  for  y 
= 0.5  agree  well  with  the  time-averaged  wave  gage 
data.  And,  the  mean  wave  profiles  for  y = 0.2 
predicted  by  both  the  DES  and  the  LES  are  seen  to  be 
close  to  the  wave  trough  visually  determined  using 
markers  and  video. 
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Figure  13.  Comparison  of  the  time-averaged  free- 
surface  elevation  predictions  by  URANS,  DES,  and 
LES  at  Fr  = 0.37  on  the  coarse  mesh  (1.8  million 
cells) 


Thus  far,  we  have  shown  and  discussed  the 
DES  predictions,  making  comparisons  with  the 
earlier  LES  predictions  whenever  possible.  Looking 
at  both  the  DES  and  the  LES  predictions,  side  by 
side,  on  the  same  fine  mesh  (4.8  million  cells)  help  us 
to  delineate  what  DES  vis-a-vis  LES  can  offer  for  the 
subject  flow.  To  further  investigate  how  URANS, 
DES,  and  LES  would  fare  with  one  another,  we 
conducted  URANS,  DES,  and  LES  computations,  on 
the  coarser  mesh  with  1.8  million  cells.  The 
resolution  of  the  coarse  mesh  is  roughly  is  half  the 
resolution  of  the  fine  mesh  used  for  the  DES 
computations  presented  so  far.  The  mean  free- 
surface  elevations  predicted  using  the  three  different 
approaches  are  shown  in  Figure  13.  Interestingly 
enough,  yet  not  surprisingly,  the  URANS  result 


obtained  with  the  realizable  k-e  model  yields  the  most 
pronounced  Kelvin  wave  system.  The  DES  captures 
the  Kelvin  waves,  although  the  wave  crests  and 
troughs  are  slightly  less  pronounced  than  the  RANS 
prediction.  The  most  remarkable  observation  is  that 
LES  on  this  coarse  mesh  fail  to  resolve  the  Kelvin 
waves. 


(c)  LES  (Kim  and  Rhee,  2006) 


Figure  14.  The  r.m.s.  wave  height  fluctuation 
predicted  by  the  present  URANS,  DES,  and  LES  at 
Fr  = 0.37  with  the  coarse  mesh 


The  capabilities  of  URANS,  DES,  and  LES 
to  reproduce  transient  aspects  of  the  free-surface 
flows  seem  to  widely  vary.  Figure  14  shows  the 
contours  of  the  r.m.s.  free-surface  elevation 
fluctuation  predicted  by  the  three  different 
approaches.  In  lieu  of  being  directly  computed,  the 
r.m.s.  wave  height  fluctuation  was  estimated  using: 


where  /"  is  the  r.m.s.  wave  height  fluctuation  and 
w'  the  r.m.s.  vertical  (z-)  velocity  fluctuation.  In  the 
earlier  study  (Kim  and  Rliee,  2006),  the  r.m.s.  wave 
height  fluctuation  estimation  using  the  LES  had  been 
found  to  be  quite  close  to  the  measured  r.m.s.  data. 
Both  the  DES  and  the  LES  predictions  are  shown  to 
yield  largely  the  same  degree  of  fluctuation.  One 
notable  difference  between  the  DES  and  the  LES 
results  is  that  the  LES  gives  a significantly  large 
r.m.s.  value  in  a larger  portion  of  the  free-surface. 
Furthermore,  the  LES  prediction  seems  to  pick  the 
fluctuation  in  the  bow  region  due  to  the  spilling 
breakers  and  the  capillary  waves.  The  DES  result  is 
also  indicative  of  the  fluctuation  at  the  bow,  yet  to  a 
less  degree.  The  URANS  gives  an  insignificant 
fluctuation  of  the  free-surface  elevation. 


Figure  15.  Locations  of  the  three  pressure  probes  in 
bow,  separation,  and  flapping  regions  marked  by 
circles 


Figures  15  and  16  show  the  time  histories  of 
local  static  pressure  signals  and  their  spectra  at  three 
selected  locations  predicted  by  the  DES  and  the  LES 
on  the  fine  mesh.  The  location  of  the  pressure  probes 
are  marked  with  a circled  “+”  symbol  in  Figure  15: 
(i)  bow-wave  region;  (ii)  separated  region;  (iii) 
shoulder  wave  region  outside  the  separated  region. 
The  spectra  shown  in  Figure  16  for  both  the  DES  and 
the  LES  indicate  that  the  pressure  signals  predicted 
are  all  largely  broad-banded,  yet  showing  distinct 
frequencies  around  which  the  power  is  concentrated. 

The  spectrum  of  the  pressure  signal  at  the 
bow  predicted  by  the  LES  has  a significant  portion  of 
power  concentrated  in  high-frequency  centered 
around  5 ~ 6 Hz  most  likely  due  to  the  weak  spilling- 
breaker  in  the  bow  waves.  The  frequency  range  is 
quite  close  to  the  measured  values  of  7 Hz  quoted  by 
Metcalf  el  a/. (2001 ) and  8.5  Hz  reported  by  Metcalf 
el  al.{ 2006)  more  recently.  The  spectra  of  the 


pressure  signal  at  the  separated  shows  pronounced 
peaks  in  the  neighborhood  of  2.0  Hz,  which  closely 
match  the  dominant  frequency  found  in  the 
experiment,  f = 1.98  Hz  (Metcalf  et  ai,  2006).  The 
flapping  region  seems  to  be  dominated  by  a low- 
frequency  mode  near /=  0.3  Hz,  as  clearly  shown  in 
both  its  spectrum.  This  frequency  is  also  close  to  the 
measured  one  (f  = 0.3  Hz).  The  DES  predictions, 
overall,  give  a lower  power  density  than  the  LES 
predictions  over  the  entire  frequency  range  for  all 
three  pressure  signals. 
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(b)  LES  (Kim  and  Rliee,  2006) 


Figure  16.  Power  spectra  of  the  local  static  pressure 
signals  recorded  at  three  different  locations  near  the 
free-surface 


Unsteady  forces  on  the  hydrofoil 


The  unsteadiness  in  the  flow  caused  by  the 
breaking  waves  and  the  large-scale  turbulent 
structures  in  different  regions  of  the  flow  leads  to 
fluctuations  in  the  resulting  forces  and  moments  on 
the  hydrofoil  surface.  The  oscillating  forces  and  the 
moments  can  adversely  affect  the  maneuvering 


characteristics  of  surface  vehicles  such  as  directional 
stability,  especially  when  the  frequency  of  the 
oscillations  in  forces  and  moments  falls  in  the 
neighborhood  of  the  natural  frequencies  of  the 
corresponding  motions  of  vehicles. 
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(b)  LES  (Kim  and  Rhee,  2006) 


Figure  17.  Time  histories  of  the  lift  and  drag 
coefficients:  (a)  DES  prediction;  (b)  LES  prediction 

Figures  17(a)  and  17(b)  show  the  time  histories  of  the 
drag  and  the  lift  coefficients  predicted  by  the  DES 
result  along  with  the  LES  prediction.  In  the  DES  and 
the  LES  results  alike,  the  drag  coefficients  oscillate 
with  a fairly  low  frequency  near /=  0.3  Hz,  which  is 
close  to  the  dominant  frequency  of  the  pressure 
signal  from  the  wave  crest  away  from  the  hydrofoil 
(the  region  of  flapping  instability)  discussed  earlier. 
The  most  remarkable  difference  between  the  DES 
and  the  LES  results  lies  with  the  lift  coefficient.  The 
lift  coefficient  predicted  by  the  LES  fluctuates  with 
much  larger  amplitude  than  that  from  the  DES.  Also 
noteworthy  in  the  history  of  the  lift  coefficients,  for 
both  the  DES  and  the  LES  alike,  is  that,  for  the 
duration  of  time  as  shown  in  the  plot,  the  time- 
averaged  lift  coefficients  predicted  by  the  DES  and 


the  LES  appreciably  deviate  from  zero.  This  aspect 
of  the  lift  coefficient  was  discussed  by  Xing  et  al. 
(2006)  and  Kim  and  Rhee  (2006). 

SUMMARY,  DISCUSSION,  AND  FUTURE 
WORK 

In  this  study,  a hybrid  URANS/LES  (DES) 
approach  was  attempted  for  the  turbulent  free-surface 
flows  past  two  surface-piercing  bodies.  Our  main 
focus  was  to  evaluate  the  DES  approach  for  turbulent 
free-surface  flows  involving  wave-breaking  and 
separation,  and  to  compare  the  results  against  LES 
and  RANS  modeling  approaches. 

• The  salient  features  of  complex  turbulent 
free-surface  flows  such  as  spilling-breaker  at 
the  bow  along  with  the  capillary  waves, 
breaking  of  the  shoulder  waves,  flow 
separation,  splashing,  and  bubble  formation 
are  largely  reproduced  by  the  DES 

• In  the  near-field,  the  DES  yields  less 
detailed  flow  features  and  a lesser  degree  of 
unsteadiness  compared  to  the  LES  results. 
This  is  primarily  due  to  the  RANS 
turbulence  model  that  is  activated  in  the 
region,  and  is  partly  due  to  different  subgrid- 
scale  modeling  in  the  DES  approach  adopted 
in  this  study. 

• The  quality  of  LES  was  found  to  degrade  as 
the  mesh  becomes  too  coarse  to  resolve 
energy  containing  eddies.  For  the  case  of 
the  hydrofoil,  under-resolved  LES  led  to  a 
failure  to  capture  the  Kelvin  waves  away 
from  the  body  where  the  mesh  stretches 
rapidly.  The  DES,  however,  performed 
commendably,  reproducing  the  Kelvin  wave 
system. 
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